In [37]:
using Plots, ComplexPhasePortrait, Interact, ApproxFun, Statistics, SpecialFunctions, LinearAlgebra
gr();
periodic_rule(n) = 2π/n*(0:(n-1)), 2π/n*ones(n)
function circle_rule(n, r)
θ = periodic_rule(n)[1]
r*exp.(im*θ), 2π*im*r/n*exp.(im*θ)
end
function ellipse_rule(n, a, b)
θ = periodic_rule(n)[1]
a*cos.(θ) + b*im*sin.(θ), 2π/n*(-a*sin.(θ) + im*b*cos.(θ))
end
Out[37]:
Dr. Sheehan Olver
s.olver@imperial.ac.uk
This lecture we cover
Quadrature rules are pairs of nodes $x_0,\ldots,x_{n-1}$ and weights $w_0,\ldots,w_{n-1}$ to approximate integrals $$ \int_a^b f(x) dx \approx \sum_{j=0}^{n-1} w_j f(x_j) $$ In this lecture we construct quadrature rules on complex contours $\gamma$ to approximate contour integrals.
The trapezium rule gives an easy approximation to integrals. On $[0,2\pi)$ for periodic $f(\theta)$, we have a simplified form:
Definition (Periodic trapezium rule) The periodic trapezium rule is the approximation
$$\int_0^{2 \pi} f(\theta) d \theta \approx {2\pi \over n} \sum_{j=0}^{n-1} f(\theta_k)$$for $\theta_j = {2 \pi j \over n}$.
The periodic trapezium rule is amazingly accurate for smooth, periodic functions:
In [2]:
f = θ -> 1/(2 + cos(θ))
errs = [((x, w) = periodic_rule(n); abs(sum(w.*f.(x)) - sum(Fun(f, 0 .. 2π)))) for n = 1:45];
plot(errs.+eps(); yscale=:log10, title="exponential convergence of n-point trapezium rule", legend=false, xlabel="n")
Out[2]:
The accuracy in integration is remarkable as the trapezoidal interpolant does not accurately approximate $f$: we can see clear differences between the functions here:
In [3]:
n=20
(x, w) = periodic_rule(n)
plot(Fun(f, 0 .. 2π); label = "integrand")
plot!(x, f.(x); label = "trapezium approximation")
Out[3]:
We can use the map that defines a contour to construct an approximation to integrals over closed contour $\gamma$:
Definition (Complex trapezium rule) The complex trapezium rule on a contour $\gamma$ (mapped from $[0,2\pi)$) is the approximation
$$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$for
$$z_j = \gamma(\theta_j) \qquad\hbox{and}\qquad w_j = {2\pi \over n} \gamma'(\theta_j)$$Example (Circle trapezium rule) On a circle $\{r e^{i \theta} : 0 \leq \theta < 2 \pi\}$, we have
$$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$for $z_j = r e^{i \theta_j}$ and $w_j = {2 \pi i r \over n} e^{i \theta_j}$.
Here we plot the quadrature points:
In [27]:
ζ, w = circle_rule(20, 1.0)
scatter(ζ; title="quadrature points", legend=false, ratio=1.0)
Out[27]:
The Circle trapezium rule is surprisingly accurate for analytic functions:
In [5]:
ζ, w = circle_rule(20, 1.0)
f = z -> cos(z)
z = 0.1+0.2im
sum(f.(ζ)./(ζ .- z).*w)/(2π*im) - f(z)
Out[5]:
In [30]:
f = z -> gamma(z)
fp = z -> gamma(z)polygamma(0,z) # exact derivative
x = 1.2
fp_fd = [(h=2.0^(-n); (f(x+h)-f(x))/h) for n = 1:50]
plot(abs.(fp_fd .- fp(x)); yscale=:log10, legend=false, title = "error of finite-difference with h=2^(-n)", xlabel="n")
Out[30]:
But the previous formula tells us that we can reduce a derivative to a contour integral. The example above shows that it's still numerically unstable, but we can deform the integration contour, to make it stable!
In [31]:
trap_fp = [((ζ, w) = circle_rule(n, 0.5);
ζ .+= x; # circle around x
sum(f.(ζ)./(ζ .- x).^2 .*w)/(2π*im)) for n=1:50]
plot(abs.(trap_fp .- fp.(x)); yscale=:log10,
title="Error of trapezium rule applied to Cauchy integral formula", xlabel="n", legend=false)
Out[31]:
We can take things further and use this to calculate the higher order derivatives, with some care taken for choosing the radius:
In [32]:
k=100
r = 1.0k
g = Fun( ζ -> exp(ζ)/(ζ - 0.1)^(k+1), Circle(0.1,r))
factorial(1.0k)/(2π*im) * sum(g) - exp(0.1)
Out[32]:
Bournemann 2011 investigates this further and optimizes the radius.
The exponential convergence of the complex trapezium rule is a consequence of $f(\gamma(t))$ being 2π-periodic, as we will see later in a few lectures:
In [33]:
θ = periodic_rule(100)[1]
plot(θ, real(f.(0.6exp.(im*θ) .+ x)./(0.5exp.(im*θ))))
Out[33]:
Example (Ellipse trapezium rule) On an ellipse $\{a \cos \theta + b i \sin \theta : 0 \leq \theta < 2 \pi\}$ we have $$\oint_\gamma f(z) dz \approx \sum_{j=0}^{n-1} w_j f(z_j)$$ for $z_j = a \cos \theta_j + b i \sin \theta_j$ and $w_j = {2 \pi \over n} (-a \sin \theta_j + i b \cos \theta_j)$.
We can use the ellipse trapezium rule in place of the circle trapzium rule and still achieve accurate results. This gives us flexibility in avoiding singularities. Consider $$ f(z) = 1/(25z^2 + 1) $$ which has poles at $\pm \I/5$. Using an ellipse allows us to design a contour that avoids these singularities:
In [34]:
scatter([1/5im,-1/5im]; label="singularities")
θ = range(0; stop=2π, length=2000)
a = 2; b= 0.1
plot!(a * cos.(θ) + im*b * sin.(θ); label="ellipse")
Out[34]:
Thus we can still use Cauchy's integral formula:
In [35]:
x = 0.1
f = z -> 1/(25z^2 + 1)
f_ellipse = [((z, w) = ellipse_rule(n, a, b); sum(f.(z)./(z.-x).*w)/(2π*im)) for n=1:1000]
plot(abs.(f_ellipse .- f(x)); yscale=:log10, title="convergence of n-point ellipse approximation", legend=false, xlabel="n")
Out[35]:
In [36]:
f = z -> sqrt(z)
function sqrtⁿ(n,z,z₀)
ret = sqrt(z₀)
c = 0.5/ret*(z-z₀)
for k=1:n
ret += c
c *= -(2k-1)/(2*(k+1)*z₀)*(z-z₀)
end
ret
end
z₀ = 0.3
n = 40
phaseplot(-2..2, -2..2, z -> sqrtⁿ.(n,z,z₀))
Out[36]:
Here we see that the approximation is vallid on the expected ellipse:
In [13]:
(ζ, w) = ellipse_rule(20, 4.0, 1.0);
ζ .= ζ .+ 4.5;
f_c = z -> sum(f.(ζ)./(ζ.-z).*w)/(2π*im)
phaseplot(-2..10, -2 .. 2, f_c)
Out[13]:
In [14]:
(ζ, w) = ellipse_rule(100, 4.0, 1.0);
ζ .= ζ .+ 4.5;
f_c = z -> sum(f.(ζ)./(ζ.-z).*w)/(2π*im)
f_c(5.3) - sqrt(5.3)
Out[14]:
Let $A \in {\mathbb C}^{d \times d}$, ${\mathbf u}_0 \in {\mathbb C}^d$ and consider the constant coefficient linear ODE
$$ {\mathbf u}'(t) = A {\mathbf u}(t)\qquad\hbox{and}\qquad {\mathbf u}(0) = {\mathbf u}_0(0) $$The solution to this is given by the matrix exponential $$ {\mathbf u}(t) = \exp(A t) {\mathbf u}_0 $$ where the matrix exponential is defined by it's Taylor series: $$ \exp(A) = \sum_{k=0}^\infty {A^k \over k!} $$ This has stability problems, so a more convenient form is as follows:
Theorem (Cauchy integral representation for matrix exponential) Let $A \in {\mathbb C}^{n \times n}$ be a diagonalizable matrix: $$ A = V \Lambda V^{-1} \qquad\hbox{for}\qquad\Lambda = \begin{pmatrix}\lambda_1 \cr &\ddots \cr && \lambda_d\end{pmatrix} $$ Let $\gamma$ be a contour that surrounds the spectrum of $A$. Then we have $$ \exp(A) = {1 \over 2 \pi i} \oint_\gamma e^z (z I - A)^{-1} dz $$
Proof Note by definition $$ \begin{align*} \exp(A) &= \sum_{k=0}^\infty {A^k \over k!} = V \sum_{k=0}^\infty {\Lambda^k \over k!}V^{-1} = V \begin{pmatrix} \E^{\lambda_1}\\ & \ddots \\ && \E^{\lambda_d} \end{pmatrix} V^{-1} \\ &= {1 \over 2\pi \I} V \begin{pmatrix} \oint_\gamma {\E^\zeta \over \zeta - \lambda_1 } \D\zeta \\ & \ddots \\ && \oint_\gamma {\E^\zeta \over \zeta - \lambda_d } \D\zeta\end{pmatrix} V^{-1} \end{align*} $$ It was important here that $\gamma$ surrounded all $\zeta_j$.
We now take out the integration from the matrix, the easiest way to see this is to apply the Trapezium rule approximation: $$ \begin{align*} V\begin{pmatrix} \oint_\gamma {\E^\zeta \over \zeta - \lambda_1 } \D\zeta \\ & \ddots \\ && \oint_\gamma {\E^\zeta \over \zeta - \lambda_d } \D\zeta\end{pmatrix} V^{-1} &= V\lim_{n \rightarrow \infty} \begin{pmatrix} \sum_{j=1}^n {w_j \E^\zeta_j \over \zeta_j - \lambda_1 } \\ & \ddots \\ && \sum_{j=1}^n {w_j \E^\zeta_j \over \zeta_j - \lambda_d } \end{pmatrix} V^{-1} \\ &= \lim_{n \rightarrow \infty} V\sum_{j=1}^n w_j \E^\zeta_j \begin{pmatrix} {1 \over \zeta_j - \lambda_1 } \\ & \ddots \\ && {1 \over \zeta_j - \lambda_d } \end{pmatrix} V^{-1} \\ &= \lim_{n \rightarrow \infty} \sum_{j=1}^n w_j \E^\zeta_j V(I \zeta_j - \Lambda)^{-1}V^{-1} \\ &= \oint_\gamma \E^\zeta V(I \zeta - \Lambda)^{-1} V^{-1} \D \zeta = \oint_\gamma \E^\zeta (I \zeta - V\Lambda V^{-1})^{-1} \D \zeta \\ & = \oint_\gamma \E^\zeta (I \zeta - A)^{-1} \D \zeta \\ \end{align*} $$
■
Demonstration we use this formula alongside the complex trapezium rule to calculate matrix exponentials. Begin by creating a random symmetric matrix (which only has real eigenvalues):
In [42]:
A = randn(5,5)
A = A + A'
λ = eigvals(A)
Out[42]:
We can now by hand create a circle that surrounds all the eigenvalues:
In [49]:
z,w = circle_rule(100,maximum(abs.(λ))+1)
plot(z)
scatter!(λ,zeros(5); label = "eigenvalues of A")
Out[49]:
Here we wrap this up into a function circle_exp that calculates the matrix exponential:
In [50]:
function circle_exp(A, n, z₀, r)
z,w = circle_rule(n,r)
z .+= z₀
ret = zero(A)
for j=1:n
ret += w[j]*exp(z[j])*inv(z[j]*I - A)
end
ret/(2π*im)
end
circle_exp(A, 100, 0, 8.0) -exp(A) |>norm
Out[50]:
In this case, it is beneficial to use an ellipse:
In [51]:
function ellipse_exp(A, n, z₀, a, b)
z,w = ellipse_rule(n,a,b)
z .+= z₀
ret = zero(A)
for j=1:n
ret += w[j]*exp(z[j])*inv(z[j]*I - A)
end
ret/(2π*im)
end
ellipse_exp(A, 50, 0, 8.0, 5.0) -exp(A) |>norm
Out[51]:
For matrices with large negative eigenvalues (For example, discretisations of the Laplacian), complex quadrature can lead to much better accuracy than Taylor series:
In [52]:
function taylor_exp(A,n)
ret = Matrix(I, size(A))
for k=1:n
ret += A^k/factorial(1.0k)
end
ret
end
B = A - 20I
taylor_exp(B, 200) -exp(B) |>norm
Out[52]:
We can use an ellpise to surround the spectrum:
In [53]:
scatter(complex.(eigvals(B)))
plot!(ellipse_rule(50,8,5)[1] .- 20)
Out[53]:
In [54]:
norm(ellipse_exp(B, 50, -20.0, 8.0, 5.0) - exp(B))
Out[54]:
If we only know $A$, how do we know how big to make the contour? Gershgorin's circle theorem gives the answer:
Theorem (Gershgorin) Let $A \in {\mathbb C}^{d \times d}$ and define $$R_k = \sum_{j=1 \atop j \neq k}^d |a_{kj}| $$ Then $$ \rho(A) \subset \bigcup_{k=1}^d \bar B(a_{kk}, R_k) $$ where $\bar B(z_0, r)$ is the closed disk of radius $r$ and $\rho(A)$ is the set of eigenvalues.
Proof
■
Demonstration Here we apply this to a particular matrix:
In [55]:
A = [1 2 3; 1 5 2; -4 1 6]
Out[55]:
The following calculates the row sums:
In [56]:
R = sum(abs.(A - Diagonal(diag(A))),dims=2)
Out[56]:
Gershgorin's theorem tells us that the spectrum lies in the union of the circles surrounding the diagonals:
In [57]:
drawcircle!(z0, R) = plot!(θ-> real(z0) + R[1]*cos(θ), θ-> imag(z0) + R[1]*sin(θ), 0, 2π, fill=(0,:red), α = 0.2, legend=false)
λ = eigvals(A)
p = plot()
for k = 1:size(A,1)
drawcircle!(A[k,k], R[k])
end
scatter!(complex.(λ); label="eigenvalues")
scatter!(complex.(diag(A)); label="diagonals")
p
Out[57]:
We can therefore use this to choose a contour big enough to surround all the circles. Here's a fairly simplistic construction for our case where everything is real:
In [58]:
z₀ = (maximum(diag(A) .+ R) + minimum(diag(A) .- R)) /2 # average edges of circle
r = max(abs.(diag(A) .- R .- z₀)..., abs.(diag(A) .+ R .- z₀)...)
plot!(Circle(z₀, r))
Out[58]:
Definition (Fourier series) On $[-\pi, \pi)$, Fourier series is an expansion of the form $$ g(\theta) = \sum_{k=-\infty}^\infty g_k e^{i k \theta} $$ where $$ g_k = {1\over 2 \pi} \int_{-\pi}^\pi g(\theta) e^{-i k \theta} d \theta $$
Definition (Laurent series) In the complex plane, _Laurent series around $z_0$ is an expansion of the form $$ f(z) = \sum{k=-\infty}^\infty f_k (z-z_0)^k $$
Lemma (Fourier series = Laurent series) On a circle in the complex plane $$ \gamma_r = \{z_0 + re^{i \theta} : -\pi \leq \theta < \pi \},$$ Laurent series of $f(z)$ around $z_0$ converges for $z \in C_r$ if the Fourier series of $g(\theta) = f(z_0 + r e^{i \theta})$ converges, and the coefficients are given by
$$ f_k = {g_k \over r^k} = {1 \over 2 \pi i} \oint_{\gamma_r} {f(\zeta) \over (\zeta - z_0)^{k+1}} d \zeta. $$Proof This follows immediately from the change of variables $z = r \E^{\I \theta} + z_0$. ■
Interestingly, analytic properties of $f$ can be used to show decaying properties in $g$:
Theorem (Decay in Fourier/Laurent coefficients) Suppose $f(z)$ is analytic in a closed annulus $A_{r,R}$ around $0$: $$A_r(z_0) = \{z : r ≤ | z| ≤ R\}$$ Then for all $k$ $$|f_k | \leq M\min\left\{{1 \over R^k} , {1 \over r^k}\right\}$$ where $M = \sup_{z \in A} |f(z)|$.
Proof This is a simple application of the ML lemma: $$ |f_k| = {1 \over 2 \pi } \left|\oint_{\gamma_1} {f(\zeta) \over \zeta^{k+1}} \D \zeta\right| = {1 \over 2 \pi}\left|\oint_{\gamma_r} {f(\zeta) \over \zeta^{k+1}} \D \zeta\right| \leq \sup_{\zeta \in \gamma_r} |f(\zeta)| R^{-k} \leq M R^{-k} $$ which similar applies deforming to $\gamma_r$. ■
Demonstration The Fourier coefficients of $$g(\theta) = {1 \over 2 - \cos \theta}$$ satisfies for $k \geq 0$ $$ |g_k| \leq {2 \over 4 - R -R^{-1}} R^{-k} $$ for all $R \leq 2 + \sqrt{3}$:
In [26]:
g =Fun(θ -> 1/(2-cos(θ)), Laurent(-π .. π))
g₊ = g.coefficients[1:2:end]
scatter(abs.(g₊); yscale=:log10, label="|g_k|", legend=:bottomleft)
R = 1.1
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")
R = 3.5
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")
R = 2+sqrt(3)-0.1
scatter!(2/(4-R-inv(R))*R.^(-(0:length(g₊))), label = "R = $R")
Out[26]: